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Abstract 

A spectral algorithm for simulating three-dimensional, incompressible, 
parallel shear flows is described. It applies to the channel, to the parallel 
boundary layer, and to other shear flows with one wall-bounded and two 
periodic directions. Representative applications to the channel and to the 
heated boundary layer are presented. 
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Nomenclature 


unit vector in x-direction 

normal wavenumber in model problem 

pressure 

time 

streamwise velocity 

normal velocity 

spanwise velocity 

streamwise coordinate 

normal coordinate 

spanwise coordinate 

scaled streamwise wavenumber 

scaled spanwise wavenumber 

pressure gradient exponent 

velocity after transforms in x and z 

velocity after transforms in x, y, and z 

specific heat 

velocity function in Falkner-Skan equations 
suction parameter 

temperature function in Falkner-Skan equations 

identity matrix 

number of grid points in x 

number of grid points in y 

number of grid points in z 

Reynolds number 



m 


w 


a 


3 

% 

5 k,£ 


S* 

K 


n 

a 

P 

v 

O 


T 

0 ) 

At 

Ay 


temperature 

Chebyshev polynomial of degree ra 
wall temperature 
streamwise wavenumber 
spanwise wavenumber 
pressure gradient parameter 
Kronecker delta function 
displacement thickness 
conductivity 

boundary layer similarity variable 

Prandtl number 

density 

kinematic viscosity 
scaled temperature 
heating parameter 
temporal frequency 
time-step 

grid spacing in normal direction 


subscripts 

k component for normal wavenumber k 

e right-hand side of linear equations 

00 free stream values 

0 equilibrium values 

& iteration parameter 


iii 




INTRODUCTION 


The development of accurate and efficient spectral methods has made 
feasible the reliable, three-dimensional simulation of the early stages of 
transition in parallel shear flows. Orszag and Kells^ pioneered the numerical 
work on channel flow. They demonstrated that linearly stable two-dimensional 
Tollmein-Schlichting waves can exhibit a strong secondary instability to 
three-dimensional disturbances of the Benney-Lin^ type for Reynolds numbers as 

Q 

low as 1000. Wray and Hussaini presented compelling evidence for the use of 
the parallel flow approximation in their Blasius boundary-layer simulation. 
Their calculation reproduced the essential features of the Kovasznay, et al.^ 
experiment up to the two-spike stage. The presence of strong secondary 
instabilities in several other linearly stable, parallel flows has been 
demonstrated by Orszag and Patera.^ An extensive comparison of numerical 
simulation with the channel flow experiments of Nishioka, et al.^ has been 
made by Kleiser and Schumann. ^ 

The calculations cited above all employed algorithms which use direct 
methods for solving the implicit equations resulting at each time-step from 
the spatial discretization. The same is true for a recent novel spectral 

O 

algorithm for curved channel flow. The cost of these direct solution methods 
is increased substantially by the addition of even minor geometric terras or 
the temporal variation of the viscosity which is essential for assessing the 
more subtle effects of heating. This is a major consideration which has led 
us to develop an algorithm which resorts to iterative methods for the solution 
of the implicit equations. 

This paper is devoted to a description of an algorithm for transition 
simulation which is based on iterative methods. It shares some common 
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features with a method developed by Morchoisne.^ The description focuses on 
the controlled boundary-layer algorithm since this problem is more involved 
than the channel and has not yet been simulated numerically. 


BASIC EQUATIONS FOR THE PARALLEL BOUNDARY LAYER 

Mean Flow 

We require the incompressible boundary-layer equations that include the 
effects of pressure gradient, suction, and/or heating. Viscous dissipation is 
neglected, and the pressure gradient and suction distributions are chosen 
compatible with similarity solutions. The Falkner-Skan equations^ are 

(pF")" + m (1 - F' 2 ) + ^ (m + 1)FF" = 0, (1) 

P 2 P 

(kH')' + o ^ ( ra + 1)FH' = 0, (2) 

1 P 

with the boundary conditions 

F(0) = F w 
F'(0) = 0 

F'(oo) = 1 (3) 

H(0) = 1 


HO) = 0 



- 3 - 


A. prime denotes differentiation with respect to the similarity variable 


n = y/u /v x . 

00 00 

The free-stream velocity is proportional to x^. 

The fluid properties are scaled with respect to their free-stream values 
(with = 1): 

m = v( T )/y 00 

K = k(T)/k 

00 

a = p C (T )/k . 

OO p oo oo 


The dimensional velocities and temperatures are related to the similarity 
variables via 

u = u F(n) 

oo 


v 


2 


J 


V u 

00 oo 

X 


[nF'(n) - (m p +l)F(n)] 


T = (T - T )H(n) + T . 

GO OO 


The parameter m p is 
parameter g by 


related to the conventional pressure 



gradient 


The wall suction is controlled by F__ and the heating effects by T and T • 

W OO 

The Reynolds number is based on free-stream velocity, viscosity, and the 
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displaceraent thickness 6 * of the uncontrolled boundary layer: 

I u » x 

6 * " i* 72 V — * 

The results in this paper pertain to water boundary layers, for which the 
appropriate empirical formulas^ are 


p (T) = 1.002 e r ^ 

r(T) = -2.303 1 1.370 + 8.36 x 10 _4 (T-293) | (T-293)/(T-164) 
k(T) = -9.901 + 0.1002 T - 1.874 x l(f 4 T 2 + 1.040 x 10 _7 T 3 

C (T) = 41.84 x [2.140 - 9.68 x 10 -3 T + 2.69 x 10 _5 T 2 - 2.42 x 10“ 8 T 3 1. (4) 

P J 

All the results below are for a free-streara temperature = 293°K. 

The numerical solutions of equations (1)— (3) were obtained by a fourth- 
order compact finite difference scheme, with a typical accuracy of 7 
significant digits. 

As noted in the introduction, the parallel flow assumption has been used 
in this work: having fixed a reference location x in the streamwise direction 
we presume that the streamwise velocity u(y) and the temperature T(y) are the 
same at all x and that the normal velocity v(y) is zero. This "mean flow" is 
not a solution to the Navier-Stokes equations. For consistency, then, we 
imagine that the Navier-Stokes equations also include a small forcing term so 
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that the "mean flow 11 is an equilibrium solution to the full nonlinear 
equations. 

In the remainder of this paper, the fluid equations will be used in 
dimensionless form, and the overbars on p" and 1c will be dropped. Velocities 
are scaled by u^, lengths by 6 , and densities by • The temperature 
variable 0 is, as usual, 

T - T 

oo 

9 = T - T * 

W “ 

The parameters which define a particular case, then, are the Reynolds number 

Re, the pressure gradient parameter 3p , the suction parameter F w , and the 

heating parameter x = T /T . 

w 00 


Linearized Equations 

Initial conditions for the numerical simulations are based on solutions 
to the Orr-Sommerf eld and Squire equations for small amplitude velocity 
perturbations. Temperature fluctuations have been ignored in these linearized 
equations, but the effect of the mean temperature upon the viscosity and mean 
flow profiles is included. Velocity perturbations are taken to be of the form 

u(x,y,z,t) - ;Cy)e l(c “ + Bz ' “ t) 

for real a and 3* The Orr— Sommerf eld and Squire eigenvalue problems are 
solved by a Chebyshev tau method.^ The Chebyshev expansions are in terms of 
the computational variable £, which is related to the physical variable y by 
the algebraic mapping 



- 6 - 


y = 


1 + 


i + g 

^e_, 

y 

max 


(5) 


The parameter y e is roughly twice the normal distance at which u = 0.5, and 
y max is the u PP er boundary in the physical domain. A typical choice for y max 
is 15, which is roughly 5 times the boundary-layer thickness. The variable E, 
has the usual Chebyshev distribution in [-1,1]: 

— cos ^ " , i=0,l,...,N . 

y y 

Roughly half the points fall within the displacement thickness, and two-thirds 
are within the boundary-layer thickness. 

Between 50 and 70 Chebyshev polynomials are used for the solution of the 
linear eigenvalue problems. The major source of error is inaccuracies in the 
numerical solution of the mean flow. Nonetheless, the eigenvalues and eigen- 
functions are reliable to 5 or 6 digits. 


Navier-Stokes Equations 

The nonlinear simulations are performed for the equations 

A 

u t + to x u = -VP + V • (pVu) + fi 
0^ + u • V0 = V • (tcV0) + g 


( 6 ) 

(7) 


V • u = 0 


( 8 ) 
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where the vorticity m = V x u, the total pressure p = p + A. |u|^, and the 
boundary conditions are 

u = 0 u = 1 

at y ** 0, at y = », (9) 

e = o 9=i 

The forcing functions are given in terms of the mean flow variables by 


f 


3 _ 

9y 


3u 0 

( p 037 _) 



"o 

(K 0 9y 


( 10 ) 


They ensure that the mean flow is a stationary solution of Equations (6) - 

( 8 ). 


NUMERICAL METHOD 

Discrete Equations 

The nonlinear three-dimensional calculations reported below are based on 
the algorithm described for two-dimensional flow in Reference 14. The spatial 
discretization is Fourier— collocation in x and z and Chebyshev collocation in 
y. The temporal discretization is backward Euler for the pressure, Crank- 
Nicolson for the normal diffusion and conduction terms, and third— or fourth- 
order Adams-Bashforth for the remaining terras in Equations (6) and (7). The 
continuity equation is enforced as a constraint at the new time level. 
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The primitive variables have two series representations which will be 
useful in this discussion. The first is 


u(x,y,z,t) 




u k ,k (y > t)e 


i(k x + k z) 

X z 


X z 


( 11 ) 


r 2 tt 
k = 7 — k , 
x L x’ 
x 


r 2 it , 

k = - — k , 
z L z 9 
z 


where 

(periodic) domain in the x and 
involves the additional series 


and L and L are the lengths of the 
X z 

z directions, respectively. The second 


°k ,k (y ’° 

x z 





m . (t) T JO, 

x* ,k z 


( 12 ) 


where % is related to y by Equation (5). Henceforth, the subscripts k x , m, 
and k z will not be written explicitly unless necessary. The collocation 
points in the periodic directions are 


X i N L ’ 
X x 


2irk 

Z k " N L » 
z z 


i = 0,1,...,N X -1 

k = 0,1, .. . jN^-l. 


(13) 


A staggered grid is employed in the normal direction. Velocities and 
temperatures are defined at the points 


? j COS N - ’ 

J y 


j - 0,1 ,...,N 


(14) 
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and the pressures at 


’j+1/2 


= cos- 


Tf( j+1/2) 


N 


0 , 1 , • • • , 1 • 


(15) 


The continuity equation is enforced at the latter points and the remaining 
equations at the former ones. The staggered grid avoids artificial pressure 
boundary conditions, precludes spurious pressure modes, and facilitates the 
solution of the algebraic equations which arise from the implicit terms in the 
time discretization. 

Chebyshev interpolation is the natural process for transferring variables 
between the grids of Equations (14) and (15). For example, consider the 
velocity component u. Let u j , for j = 0,1,..., Ny, denote its values at the 
points (14). The Chebyshev coefficients are given by the usual quadrature 
rule 


N 

« ■ — *rYz A -\ t ) — — yV x u. 

N c “ J J J N c “ J J 


IMTJ i 

cos , m » 0,1,... ,N 


y m j=0 


y m j=0 


y 

(16) 


where 


c = 
m 


2 m = 0 or N 


1 1 < m < N 


The interpolated values of u are 


N -1 

y 


u 


j+1/2 ’ X "m V 5 j+l/2 ) - X "» C ° S I±L 1 ' N ' >1,m • J ■ °> 1 V 1 ' U7) 

m=0 m=0 ^ 


N -1 

y 


T (5 4j . 1/9 ) 
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[Note Chat T N (5j +1 / 2 ^ = 0 for j* 0 * 1 * • • • » N y _1 * ] The Fast Fourier Transform 
(FFT) may be used to evaluate both sums (16) and (17). The less familiar sura 
in Equation (17) over the odd cosines may be handled by the technique 
presented in Appendix C of Reference 15. 

The temporal discretization of Equations ( 6 )— ( 8 ) leads, after a Fourier 
transform in x and z, to an implicit system of the form 


"n+1 

u 


, "n+1. a. M An+1 

" ( bU y >y + lk x ™ 


u 

e 


~n+l ‘ ^n+L , A n+1 A „ r 

v - (bv ) + Yq = v at g. 

y y y e j 


n+1 n+lv , . <> n+1 

w - (bw ) + lk Yq = w 

y y z e 


(18) 


, i * n+1 jl n+1 . , jl n+1 n r / 1Q \ 

-ik y* u - - lk y*w =0 at £, , - /0 (19) 

x y z j+1/2 

along with Dirichlet boundary conditions on the velocities. An * denotes the 
complex conjugate. The coefficient b = ^ At ^avg( y»t)» where the last term is 
the average value of y n (x,y,z,t) at fixed y and t. The pressure has been 

A ^ A 

included in terras of the scaled variable q = (— )P, where y is a complex 

constant whose role is explained below. Similarly, the temperature equation 

is 

e n+1 - (d0 n+1 ) =0 at g. (20) 

y y e j 

subject to Dirichlet boundary conditions, where d = y At K^ y g(y,t). The 

right-hand-sides of these equations contain the explicit terras in the temporal 
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discretization. Let 


tt - /* n+ l 

U = (u Q , 


'n+1 


u. 


» • • • > 


A n+1 N 

u m )» 


Q 


, A n+l "n+1 *n+l 

W l/2’ q 3/2 * * * * * q N -1/2' 

y 


and define V, W, and 0 similarly to U. Let and be diagonal matrices 
with the elements of b and d on their respective diagonals. Let the effect of 
Equations (16) and (17) on U be denoted by A+. and the reverse interpolation 
procedure (for Q) by Aq. Finally, let M denote the matrix which represents 
Chebyshev differentiation in the y direction. Then Equations (18) to (20) 
reduce to the algebraic set 


(I - MD,M)U + ik yA^Q = U 
b x 0 e 

(I - MD M)V + yMA^Q = V (21) 

b (J e v 7 

(I - MD M)W + ik yA^Q = W , 
b z 0 e 


-ik y*A U - y*A.MV - ik y*A W =0 (22) 

XT T 2 t 

(I - MD^M)0 = 0^ (23) 

where the first and last rows of (21) and (23) are replaced by the boundary 
conditions. Clearly, the equations for each pair (& x ,k z ) are independent. 
Moreover, Equation (23) is not coupled to Equations (21) and (22) and thus may 
be solved separately. The matrices Aq, A + , and M are full. Except in special 
cases the direct solution of these equations is not practical. 
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The iterative solution of Equation (23) is straightforward: simply 
precondition the system by a finite-difference approximation on the Chebyshev 
grid^® and apply a standard iterative method^ such as Richardson, 1® Chebyshev 
acceleration, minimum residual,!® or even multigrid* ^ The finite- 
difference system is tridiagonal and positive definite* 

The key to this algorithm is the solution of the system (21)-(22). A 
simplified model problem, discussed in the following sub-section, is 
instructive. 


Model Problem Discussion 

Suppose that the boundary conditions in the normal direction are periodic 
instead of Dirichlet and that the viscosity, i.e., b, is constant* Replace 
the Chebyshev discretization with a Fourier one, on, say [0,2 tt]* Then the 
vertical collocation points are 



j-0 , 1 , . • • ,Ny 1 


_ 2tt( j+1/2) 
y j+l/2 N 

J y 


j=o,i 


f 


N -1. 

y 


The fully discrete equations may be cast in a form analogous to (21)-(22), 
where now 


M - C*DC 0 


MD b M - bCgD C Q 


A 0 = c o c + 


( 24 ) 
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A + ■ c i c 0 


and 


<c o> M 



(c A,j 



iky 

e 


j+1/2 

> 


j= 0,1,...,N -1 

k= -N y /2, -N y /2+l,...,N /2-1 




k,A = -N /2, -N /2+l,...,N /2-1. 

y y y 


Thus, we have for the spectral equations 


(I - bC* D 2 C Q )U + ik x yC* C + Q = U 0 

(I - bC* D 2 C Q )V + Y C* DC + Q = V e (25) 

(I - bC* D 2 C 0 )W + ik z Y C* C + Q = W e 


-ik x y*C* C Q U - y*C* D*C q V - ik z y*C* C Q W = 0. 
This can be written as the system 

LX = B 


(26) 


( 27 ) 


where, for instance, X = (U, V, W, Q). Now let U, = C„ U. 0 = C 0 

k 0 * x k + x 

X = RX, and L = RLR* where 
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( C 0 0 o o 

0 c 0 0 0 

0 0 C Q 0 

0 0 0 c + 

After a permutation of the rows and columns of L, we obtain a block diagonal 
matrix with blocks 



L k \ = B k 


(28) 


/v . /v 'N. rw v 

\ = K> v w k> 


L k = 


CM 

+ 
1 — H 

0 

0 


k 

° 

1+bk 2 

0 

iky 

A 

I (29) 

l 0 

0 

l+bk 

ik y 

/ 

\ A 


A 

Z i 

! 

\-i V * 

-iky* 

-ik y* 
z 

0 t 

f 

• 


Consider now a finite-difference approximation to this model problem. Let E 
denote the foward shift operator subject to periodic boundary conditions. 
Then Equations (21)-(22) become 


[i ^-y (E - 21 + E _1 ) ]u + ik Y y (E + I)Q = u 

(Ay) 2 X 6 


[I ^~2 (E - 21 + E L]V + X. ( E - I)Q = V e 

(Ay) y 


[I (E - 21 + E _1 ) ]w + ikY 1 (E + I)Q = W , 

(Ay) 2 


- ik y* i (I + E -1 )U - -7— (I - E 1 )V - ik y* y (I + E _1 )W 
x 2 Ay z 2 


0 . 


(30) 


(31) 
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Denote the matrix which represents the left-hand-side by H. This system can 
be reduced to block-diagonal form by the same transformation that was used for 
the spectral operator. The result can be written 


/V 

w - \ 


/l+bk 2 s 2 

0 

0 

03 

>- 

1 

•H 

l 0 

2 2 

1+bk s 

0 

ikys \ 

\ -° 

0 

2 2 

1+bk s 

A 

ik zY a 1 

\-ik 

X 

-iky*s 

-ik Y*a 
z 

0 / 


where 

= sin(kAy/2) 
(kAy/2) ’ 


a = cos(-^^) . 


(32) 


(33) 


(34) 


(35) 


The relevant range is |kAy| < ir. 

If a and s were identically one, then the preconditioning would be 

perfect. In any case, the derivative terms cause no serious problem 

for (2 /tt) < s < 1. It is the averaging operator a which is a source of 

potential difficulty. As |kAy| + tt, the averaging becomes useless. We 

anticipate difficulty only in circumstances for which k and/or k are large 

x z 

relative to the reciprocal of the grid spacing in y. 

A preconditioned Richardson iteration for Equation (27) reads 


X «- X - Q H -1 (B - LX) 


(36) 
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where ft is an iteration parameter. The choice of ft and indeed the convergence 
properties of this scheme depend upon the eigenvalues of H“^L. A (complex) 

Chebyshev acceleration of the basic Richardson scheme can be devised provided 

-1 20 
that all the eigenvalues of H L L have positive real parts. 

The eigenvalues of the model problem are especially easy to obtain: 

since H and L were block-diagonalized by the same transformation, we need only 

compute the eigenvalues of h"*L, for k=l ,2, . . . ,N /2. Some results are shown 

1C K. J 

A 

in Figure 1. In these calculations k z has been set to zero. Similar 

A A 

calculations for both k and k non-zero lead to qualitatively similar 
results. Figure 1(a) portrays the easiest of these 6 cases for the iterative 
scheme. Most of the eigenvalues are near unity, and they are located near the 

A 

real axis between 1 and tt/ 2. The eigenvalue spread in part (b) for k^ = 30 is 
much larger. Nevertheless, the real parts of the eigenvalues are safely 

greater than zero. A comparison of parts (b), (c), and (d) reveals that for 

A 

fixed k the eigenvalue spread is reduced as the vertical resolution is 
increased. The heuristic explanation for this welcome behavior is that as 
increases, the eigenvalues of the first derivative operator become more 
important than those of the averaging operator. In actual numerical 

simulations the number of points in the x and y directions is likely to be 
nearly the same. Part (d) corresponds to the worst case that would arise in a 
64^ calculation. Parts (e) and (f) show the eigenvalue spectrum for a 
situation in which the viscosity is quite considerable. The major difference 
from the previous cases is the presence of additional eigenvalues along the 
real axis as large as tt 2 /4. This is characteristic of preconditionings of 
second-derivative spectral operators. 
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Clearly, the major shortcoming of this particular preconditioning is its 
treatment of the averaging operator. If the 2/3 rule is used to de-alias the 
collocation approximation in the normal direction, then the averaging operator 
is well-behaved. In this event, |kAy| < 2tt/3, so that 1/2 < a < 1. For the 
case shown in Figure 1(b), the largest 6 eigenvalues disappear. With the use 
of the 2/3 rule, this case becomes quite manageable. 

In the next sub-section, we present numerical evidence that the model 
problem predicts very well the eigenvalues of the system of real interest. 
Already the model problem suggests that the major shortcoming of the 
preconditioning is its treatment of the averaging operator. Fortunately, the 
model problem provides a ready tool with which to check the effectiveness of 
alternative preconditionings. 

A Minimum Residual Iterative Scheme 

The actual system that must be solved is given by Equations (21) and 
(22). It clearly can be written in the form of Equation (27) by an obvious 
adaptation of the notation of the previous sub-section. Likewise, let H 
represent the finite difference counterpart of L on the Chebyshev staggered 
grid. 

The eigenvalues of some channel flow cases are shown in Figure 2. There 
are no significant differences between these results and those for the model 
problem. The eigenvalue distribution for the boundary layer is slightly 
broader than for the channel. See Reference 14 for some boundary layer 
eigenvalues. 

A useful alternative to the Chebyshev iterative method discussed in the 
preceding sub-section is a class of variational iterative methods. 18 The 
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simplest such scheme is called the minimum residual (MR) method. The 
preconditioned version begins with an initial guess X Q> the initial 
residual R Q = B - LX q and the initial direction Q 0 , determined by HQ 0 = R 0 > 
and proceeds according to 

- ( V l V /<L V L V 



Like the Richardson scheme, this method requires one evaluation of the 
spectral operator (for LQ n ) and one solution of the implicit finite-difference 
system (for Q n +j) per iteration. A sufficient condition for convergence is 
that the symmetric part of LH~^ be positive definite. The constant y that 
appears in Equations (21) and (22) is used to ensure that this condition is 
met. (One can easily show that y has no effect upon the eigenvalues of H“^L, 
but that it does influence the symmetric part of LH"~*.) 

Let us return momentarily to the model problem. For y = 1, the extreme 
eigenvalues of the symmetric part of LH*” 1 are 1.59 and 0.71 for the case shown 
in Figure 1(a) whereas they are 18.8 and -15.6 for Figure 1(b). The MR method 

i 2 *2 

will clearly fail for the latter case. However, for y = 1/ /k + k these 

V x z 

latter eigenvalues improve to 2.55 and 0.43. The importance of y is 
evident. It is even more important when there is appreciable diffusion. 
The y=l extreme eigenvalues are 474 and -471 for Figure 1(f). The choice 
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Y = [l + dk ]/ 

1 max J 


k 2 + k 2 + k 2 , 

x z max 


where k max = N y /2 leads to 2,98 and 0.98. This scaling is based on balancing 
the norms of the diffusion and gradient operators. The principle applies to 
the actual channel and boundary layer problems as well. 

There is clearly the prospect of future improvements in the 
preconditioning and scaling. One intriguing scheme is suggested by the 
observation that if a non-staggered grid were used for the pressure, then the 
preconditioning problems would shift from the averaging operator (which would 
become the identity) to the first derivative operator. Perhaps one should 
employ a scheme which alternates iterations on a staggered grid with ones on a 
non-staggered grid. 

The preconditioning matrix H is block-tridiagonal. Note that Equation 
(27) can be separated into 2 independent real systems. Although we have no 
proof that the linear system HQ n+1 = from Equation (37) can be inverted 
without pivoting, we have yet to encounter a case which requires it for 
numerical stability, provided that the equations are ordered as suggested in 
Reference 14 and scaled as suggested above. We have, of course, made 
comparison with calculations performed with and without pivoting. We have 
also made several production runs in both 32-bit and 64-bit arithmetic and 
found 5 digit agreement. Furthermore, substantial round-off errors arising 
during the solution of the linear equations should prevent the iterative 
scheme from converging and this has not been observed. (The MR method, when 
convergent, has the property that the residual can be reduced to an 
arbitrarily small level, regardless of the precision of the machine.) 
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Time Discretization 

This algorithm does not resort to splitting. * Although a splitting 
scheme would permit a simpler iterative method, this approach appears 
inadvisable for a general code in view of the demonstration by Marcus^ 1 for 
Taylor-Couette flow that the standard splitting scheme produces errors that do 
not vanish as At tends to zero. 

The Crank-Nicolson treatment of the mean vertical diffusion term is 
standard and is essential for practical calculations at low Reynolds number. 
If needed, a serai-implicit treatment of horizontal diffusion can be readily 
incorporated into the algorithm. Moreover, the mean streamwise advection may 
also be treated semi-implicitly • 

A backwards Euler treatment of the pressure appears to be all that is 
warranted. This variable merely serves as a constraint for enforcing the 
incompressibility condition. The pressure term has the character of an 
advection term with an infinite speed. The potential hazards of Crank- 
Nicolson for advection terms are well-known. We ourselves have encountered 
some difficulty with the use of Crank-Nicolson on the pressure in problems 
characterized by rapid decay of the interesting components of the solution. 
One such case is the Stokes layer in channel flow with one oscillating wall. 
The computed velocity decays properly but the pressure does not. It quickly 
attains an amplitude which does not vary with time and it changes sign every 
time step. No such difficulty arises when the backward Euler scheme is used 
for the pressure. Moreover, the velocity agrees with that of the former 
calculation. We have made numerous comparisons of the two pressure treatments 
in more conventional problems in which the solution grows or only slowly 
decays. In no case has there been any detectable difference in the 
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velocities* Evidently, the overall time accuracy of the algorithm is not 
degraded by the backward Euler pressure scheme* 

Many transition and turbulence simulations have used second-order Adams- 
Bashforth (AB2) for the advection term. Both this method and explicit second- 
order Runge-Kutta (RK2) methods for a Fourier spatial discretization of 
advection suffer from weak instability* For a given spatial grid the fully 
discrete equations have a parasitic solution with a positive growth rate which 
is proportional to At (in terms of the physical time) • This means that the 
computed solution is destined to blow up if integrated long enough* The 
useable time interval can be increased by reducing At, but this can be 
burdensome for long time integrations* 

Higher-order time-stepping methods have the advantage of asymptotic 
stability as well as improved accuracy* A practical advection stability 
condition has the form 


At < c 


-**_ + + -JlL_ 

u W Jmax w 

max max 


(38) 


where c is a CFL parameter. The formal limits on c are 0.23 for AB3, 0.13 for 
AB4, 0.55 for RK3, and 1.27 for RK4. In terms of execution time for a 
calculation operated at the stability limit, AB3 is favored over RK3 and RK4 
over AB4. 

The Navier-Stokes algorithm with a conventional third-order advection 
scheme requires storage for 10 variables. (The pressure is needed only at the 
latest time level.) Williamson^ has catalogued numerous low storage Runge- 
Kutta schemes that, when modified to the present application with its 
additional implicit terms, require storage for only 7 variables. For the 
equation 
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u = F + G, 
t 


where the terms in F are treated by RK3 and those in G by Crank-Nicolson, 
such a scheme reads 

H. = AtF 
1 o 


u, = u + I H. + 1 At(G + G ) 
1 o 3 1 6 o 1 


H 


u 


H 


u 


2 


2 


3 


3 


AtF x 


U 1 + 16 H 2 + 24 At ^ G l + G 2 ^ 


. t, 153 „ 
AtF 2 “ 128 H 2 


u 2 + T5 H 3 + ¥ At(G 2 + G 3> 


(39) 


where u Q = u n and U 3 = u n+ ^. The inclusion of G makes this scheme formally 
second-order accurate for advection. However, in our own calculations for low 
viscosity flows we have found it to have errors at most 50% greater than those 
for a true third-order advection scheme. Moreover, the errors decrease by a 
factor of 8 when At is halved, even down to an accuracy level of 6 significant 
digits. The scheme (39) is to be preferred when storage is at a premium or 
the I/O costs are substantial. 
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Implementation 

Both the channel and boundary-layer versions of this algorithm have been 
implemented on a Control Data Corporation Cyber 175. The channel code has 
recently been made operational on the Control Data Corporation VPS 32 (a 
vector processor which is architecturally similar to a two-pipe Cyber 205, but 
expanded to over 16 million 64-bit words of memory). 

No special coding was used for the scalar Cyber 175 code, except for the 
assembly language FFT's. The code has been used for calculations on 
collocation grids as large as 16 x 32 x 8. A typical channel application 
requires 3 msec per time step per grid point for a convergence criterion on 
the iterative scheme sufficient to ensure that the velocity field is 
divergence free to better than 1 part in 10 10 . (If the 2/3 rule is used to 
de-alias the horizontal directions, then this time is reduced substantially—— 
roughly by half — since most of the CPU time is spent in solving the implicit 
equations.) The run time does depend on the amplitude of the disturbances. 
It varies by perhaps a factor of 2 in either direction from the figure cited 
above. The boundary-layer code, with fixed temperature, takes perhaps 20% 
more time per step. 

The VPS 32 code has been implemented entirely in Cyber 200 Fortran. The 
FFT s were written by the authors using the guidelines given by Temperton.^ 
They currently incorporate radix 2 and radix 3. Vectorization of the implicit 
equation solution was achieved by solving for many pairs (k x , k z ) 
simultaneously with no pivoting. Typically 1/4 of all the pairs were solved 

o 

for at a time. On a 32 grid the vector lengths for the block-tridiagonal 
inversions and for many parts of the residual calculation are only 136. The 
VPS 32 requires somewhat longer vectors . to operate near its peak capacity. 
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For a 64 3 grid, however, the vector length is over 500, which is much more 
acceptable. A typical calculation on a 32 3 grid requires 2.5 sec per step, 
whereas the corresponding 64 3 calculations take 10 sec per step. This 
increase in time by a factor of 4 rather than 8 is indicative of the improved 
performance of the VPS 32 as the vector length increases. (With aliasing 
control these times are halved.) 

We have also implemented a version of this code which uses finite 
differences in the normal direction. It takes 0.5 sec per step on a 32 3 grid 
and 3 sec per step on a 64 3 grid. These calculations proceed at a sustained 
speed of over 90 MFLOPS in 64-bit arithmetic. Clearly, the CPU time penalty 
for spectral resolution in the normal direction is not unduly severe. 


APPLICATIONS 

Channel Flow 

The usual scaling for the channel is employed. Lengths are scaled by 
channel half-width and velocities by the centerline velocity of the mean 
flow. A uniform density of unity is presumed. The Reynolds number is based 
on channel half-width and the mean centerline velocity. 

Parallel shear flows admit one set of linear waves which are solutions to 
the Orr-Somraerf eld equation and another set which are solutions to the 
unforced Squire equation. The linear waves for the channel can be further 
delineated into wall modes and center modes. The former have phase speeds of 
roughly 0.3, and their most significant details are located near the walls. 
The Orr-Sommerf eld wall modes are, of course, the familiar Tollmein- 
Schlichting (TS) waves. Although the Squire wall modes (SW) are themselves 
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linearly stable, weakly nonlinear theories have indicated that they may 
interact with TS waves, even some linearly stable ones, in such a way as to 
trigger instability* The center inodes, both the Orr-Sommerf eld ones (OSC) and 
the Squires ones (SC), have phase speeds roughly 0.9, and their energy is 
concentrated near the center of the channel. 

Because of the availability of experimental^ and computational^ results, 
an obvious test case is ‘ Re « 5000. Table I summarizes the properties of 
several linear modes. Columns 4 and 5 list the real and imaginary parts of 
the temporal eigenvalues m, as determined by the numerical solution of the 
linear eigenvalue problem. The last column in Table I lists the growth rates 
estimated from numerical simulations using as initial conditions just the mean 
flow plus a linear mode at 0.001% amplitude, where the eigenfunctions are 
normalized so that they have a maximum strearawise velocity component of 1. 
The calculated growth rates were obtained from a least squares fit to the 
perturbation kinetic energy over the time interval [0,50]. These particular 
calculations used = 48 and the AB3 scheme with At = 0.100. This seems to 
be adequate evidence for the consistency of the present semi-implicit, 
unsplit, collocation, staggered grid spectral algorithm. 

The time history of the kinetic energy in selected harmonics is displayed 
in Figure 3 for a more interesting calculation. This Re = 5000 run on a 
16 x 32 x 16 grid began with the two-dimensional TS wave of Table I at 5% 

amplitude and the usual Benney-Lin combination of two oblique TS waves at a 

combined amplitude of 0.1%. The lowest harmonics retained in the calculations 
were 1.12 in x and 2.00 in z. The time step was At = 0.05. The individual 
harmonics are labeled by integers which denote the wavenumbers relative to 

those of the three-dimensional TS waves. The modal energies are measured 
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relative to the energy in the mean flow and include contributions from 
symmetrical components; e.g., the (1,-1), (-1,1), and (-1,-1) components are 

included in the harmonic labelled (1,1). Nonlinear effects produce the 
familiar secondary instability which is indicated by the growth of the (1,1) 
mode. 

Heated, Parallel Boundary Layer 

The controlled, incompressible boundary layer has not yet been explored 

by three-dimensional numerical simulation. The principal LFC techniques can 

be incorporated into the parallel boundary layer approximation as described in 

the second section. We report here some of the tests of the present algorithm 

which were performed with the low resolution Cyber 175 code with temporally 

frozen temperature. 

£ 

The case Re = 8950 has been singled out in the past for investigations 
based on linear theory.^ The strongest instability occurs for a = 0.158266, 

3 = 0 and has a complex frequency = 0.036797 + i 0.003550. Each of the 
principal LFC techniques can completely stabilize flow at this Reynolds 
number, at least linearly. Table II lists the amount of control required to 
reduce the growth rate of the strongest instability to roughly 0.0001. This 
table includes the wavenumbers and frequencies of the least stable (two- 
dimensional) waves, as well as the frequency of the three-dimensional wave 
with 3 = ot. The last column is indicative of the accuracy of the numerical 
simulations. It lists the growth rates estimated by a least squares fit to 


*Bushnell, D. M., M. Y. Hussaini, and T. A. Zang, "Sensitivity of LFC 
Techniques in the Nonlinear Regime," to appear. 
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the perturbation kinetic energy for a calculation starting with an initial TS 
wave of 0.01% amplitude and lasting until t = 500. The time-step was At = 1 
and the vertical resolution was N z = 40. The integration time covers well 
over two periods of each wave. The discrepancy between the growth rates 
estimated from the simulation and those predicted by linear theory is of the 
same order as the accuracy of the latter. 

Figure 4 demonstrates that the heated boundary layer is susceptible to 
the secondary instability. This calculation was performed on a 16 x 32 x 8 
grid with a time step At = 0.5. The usual initial three-dimensional Benney— 
Lin mode had an amplitude of 0.01% and the two-dimensional TS wave started at 
the 5% level. The time histories of the harmonic components are very similar 
to those for the secondary instability of the channel wall modes. 

Numerous other examples for both the channel and the boundary layer may 
be found in reference 24. 


CONCLUDING REMARKS 

This paper has been devoted to a detailed description of a fully spectral 
algorithm that we are presently employing to investigate stability and 
transition in controlled, parallel boundary layers. This method is 
particularly well-suited to accounting for the subtle, time-dependent effects 
of heating and of temperature fluctuations upon the viscosity and 
conductivity. 

This algorithm may be extended to other applications as well. It can, 
for example, be used for large-eddy simulations or for other simulations of 
turbulence which employ a spatially and temporally varying eddy viscosity. 
Moreover, it may also be applied to compressible flow problems. 
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Table I. 

Some 

Channel Modes for 

Re = 5000 


Mode 

a 

e 

03 

r 

“i 

calc 

TS 2-D 

1.12 

0.00 

0.315563 

-0.002783 

-0.002783 

IS 3-D 

1.12 

2.00 

0.364473 

-0.076227 

-0.076234 

SW 3-D 

0.56 

2.00 

0.125129 

-0.069908 

-0.069889 

OSC 2-D 

1.12 

0.00 

1.067058 

-0.052467 

-0.052466 

OSC 3-D 

1.12 

2.00 

1.066915 

-0.051005 

-0.050997 


Table 

II 

• Some 

Controlled Boundary-Layer Modes 

for Re = 

8950 

Control 

Mode 

a 

6 

03 

r 

03. 

l 

W i | calc 

3 = 0.55 

T) 

TS 

2-D 

0.167675 

0.00 

0.037384 

0.000095 

0.000096 

r 

TS 

3-D 

0.167675 

0.167675 

0.040948 

-0.001012 

-0.001028 

F = 0.895 
w 

TS 

2-D 

0.162057 

0.00 

0.036207 

0.000093 

0.000093 

TS 

3-D 

0.162057 

0.162057 

0.039742 

-0.000968 

-0.000993 

T = 1.10 

TS 

2-D 

0.149937 

0.00 

0.029337 

0.000093 

0.000097 


TS 

3-D 

0.149937 

0.149937 

0.032106 

-0.000798 

-0.000793 
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Eigenvalues^ of H~^ L for channel flow, 
values of k and k . 
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Figure 3. Harmonic history for a Re = 5000 channel flow simulation starting with a 
5% two-dimensional TS wave and a 0.1% combination of 2 three-dimensional 
TS waves. The harmonics are labelled by their respective streamwise and 
spanwise wavenumbers relative to those of the three-dimensional TS wave 
with a = 1.12 and 6 = 2.00. 
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Figure 4. Harmonic history for a Re = 8950 heated boundary layer simulation 
starting with a 5% two-dimensional TS wave and a 0.01% combination of 2 
three-dimensional TS waves. The harmonics are labelled by their 
respective streamwise and spanwise wavenumbers relative to those of the 
three-dimensional TS wave with a = 0.150 and S - 0.150. 
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